home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / examples / doc / sigprc14 < prev    next >
Text File  |  1997-07-08  |  1KB  |  46 lines
  1. ; This batch file creates a plot of the impulse response and the 
  2. ; frequency response of a notch filter, as described in the section
  3. ; on Infinite Impulse Response filters in Chapter 13, "Signal Processing",
  4. ; of _Using IDL_.
  5.  
  6. @sigprc13.bat   ; Load the coefficients for the notch filter.
  7.  
  8. na = N_ELEMENTS(A)-1 ; degree of denominator polynomial
  9. nb = N_ELEMENTS(B)-1 ; degree of numerator polynomial
  10. N = 1024L
  11.  
  12. ; set the input u to an impulse
  13.  
  14. U = FLTARR(N)
  15. U(0) = FLOAT(N)
  16.  
  17. Y = FLTARR(N)
  18. Y(0) = B(2)*U(0) / A(na)
  19.  
  20. FOR K = 1, N-1 DO $ ; compute filtered signal recursively
  21.    Y(K) = ( TOTAL( B(nb-K>0:nb  )*U(K-nb>0:K  ) ) $
  22.           - TOTAL( A(na-K>0:na-1)*Y(K-na>0:K-1) ) ) / A(na)
  23.  
  24. V = FFT(Y) ; compute spectrum v
  25.  
  26. F = FINDGEN(N/2+1) / (N*delt) ; f = [0.0, 1.0/(N*delt), ... , 1.0/(2.0*delt)]
  27. mag = ABS(V(0:N/2)); magnitude of first half of v
  28. phi = ATAN(V(0:N/2)) ; phase of first half of v
  29.  
  30. ; log plots of magnitude in dB and phase in degrees
  31.  
  32. !P.MULTI = [0, 1, 2, 0, 0] ; set up for two plots in window
  33.  
  34. PLOT, F, 20*ALOG10(mag), YTITLE='Magnitude in dB', $
  35.     XTITLE='Frequency in cycles / second', $
  36.     /XLOG, XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1, $
  37.     TITLE='Frequency Response Function of b(z)/a(z)'
  38.  
  39. PLOT, F, phi/!DTOR, YTITLE='Phase in degrees', $
  40.     YRANGE=[-180,180], YSTYLE=1, YTICKS=4, YMINOR=3, $
  41.     XTITLE='Frequency in cycles / second', /XLOG, $
  42.     XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1
  43.  
  44. !P.MULTI = 0
  45.  
  46.